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Abstract— A simple numerical procedure for estimating the stochastic 
robustness of a linear time-invariant system is described. Monte Carlo 
evaluation of the system's eigenvalues allows the probability of instabil- 
ity and (he related stochastic root locus to be estimated. This analysis 
approach treats not only Gaussian parameter uncertainties but non- 
Gaussian cases, including uncertain-but-bounded variations. Confidence 
intervals for the scalar probability of instability address computational 
issues inherent in Monte Carlo simulation. Trivial extensions of the 
procedure admit consideration of alternate discriminants; thus, the 
probabilities that stipulated degrees of instability will be exceeded or 
that closed-loop roots will leave desirable regions can also be estimated. 
Results are particularly amenable to graphical presentation. 

I. Introduction 

Control system robustness is defined as the ability to maintain 
satisfactory stability or performance characteristics in the presence 
of all conceivable system parameter variations. While assured ro- 
bustness may be viewed as an alternative to gain adaptation or 
scheduling to accommodate known parameter variations, more often 
it is seen as protection against uncertainties in plant specification. 
Consequently, a statistical description of control system robustness 
is consistent with what may be known about the structure and 
parameters of the plant's dynamic model. 

Guaranteeing robustness has long been a design objective of 
control system analysis, although in most instances, insensitivity to 
parameter variations has been treated as a deterministic problem 
(see [1] for a comprehensive presentation of both classical and 
modem robust control). Stability (gain and phase) margins are 
useful concepts for designing robust single-input /single-output sys- 
tems, addressing disturbance rejection and other performance gois, 
and they are amenable to the manual graphical procedures that 
preceded the widespread use of computers. With the help of com- 
puters, singular- value analysis has extended the frequency -domain 
approach to multiinput/multiourput systems (e.g., [2], [3]); how- 
ever, guaranteed stability -bound estimates often are unduly conser- 
vative, and the relationship to parameter variations in the physical 
system is weak. Structured singular-value analysis [4] reduces this 
conservatism somewhat, and alternate treatments of structured pa- 
rameter variations have been proposed (e.g., [5] -[7]), although 
these approaches remain deterministic. Reference [8] uses the term 
“stochastic robustness" to describe a stability bound based on 
Lyapunov methods and parameter perturbations that are modeled as 
stochastic sequences. This is a deterministic stability bound ex- 
pressed in terms of the norm of a vector of noise variances. 
Elements of “stochastic stability" [9] have application to robustness 
but have yet to be presented in that context. 

The notion of probability of instability , which is central to the 
analysis of stochastic robustness, was introduced in [10], with 
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application to the robustness of the Space Shuttle's flight control 
system, and it is further described in [11]- [14]. This method 
determines the stochastic robustness of a linear time-invariant 
system by the probability distributions of closed-loop eigenvalues, 
given the statistics of the variable parameters in the plant's dynamic 
model. The probability that all of these eigenvalues lie in the open 
left-half 5 plane is the scalar measure of robustness. 

With the advent of graphics workstations, the stochastic robust- 
ness of a system is easily computed by Monte Carlo simulation, and 
results can be displayed pictorially, providing insight into otherwise 
hidden robustness properties of the system. The method is computa- 
tionally simple, requiring only matrix manipulation and eigenvalue 
computation, and it is inherently nonconservative, given a large 
enough sample space. Furthermore, the analysis of stochastic ro- 
bustness is a logical adjunct to parameter-space control design 
methods [15]- [18]. Details of the approach and an example are 
given in the sequel. 

n. Probability of Instability 

Consider a linear time-invariant (LTI) system subject to LTI 
control 

*(0 = r(p)x(t) + G(p)u(t) (i) 

y{t) = H(p)x{t) ( 2 ) 

u (0 = "<•(') ~ CH(p)x{t) ( 3 ) 

x(t), u(t ), y(t), and p are state, control, output, and parameter 
vectors of dimension n , m , q , and r, respectively, accompanied by 
conformable dynamic, control, and output matrices F, G, and H , 
which may be arbitrary functions of p. u c (t) is a command input 
vector, and, for simplicity, the (m x n) control gain matrix C is 
assumed to be known without error. The n eigenvalues, X, = o i + 
y'o),, / = 1 to n, of the matrix [ F(p ) - G( p)CH( p)] determine 
closed-loop stability and can be determined as the roots of the 
determinant equation 

UL - [T(p) - G(p)CH(p)] I = 0 (4) 

where s is a complex operator and l n is the (n x n) identity 
matrix. While the explicit relationship between parameters and 
eigenvalues is complicated, estimating the probability of instability 
of the closed-loop system from repeated eigenvalue calculation is a 
straightforward task. Putting aside the mathematical intricacies, note 
that the probability of stability plus the probability of instability is 
one. Since stability requires all the roots to be in the open left-half s 
plane, while instability results from even a single right-half s plane 
root, we may write 

/ o 

pr {a) da (5) 

- OD 

where a is an n- vector of the real parts of the system's eigenvalues 
(X = o + yw), pr(a) is the joint probability density function of a 
(unknown analytically), and the integral that defines the probability 
of stability is evaluated over the space of individual components of 
a. Denoting the probability density function of p aLS pr(p), (4) is 
evaluated J times with each element of p r j — 1 to 7, specified by 
a random-number generator whose individual outputs are shaped by 
pr(p). This Monte Carlo evaluation of the probability of stability 
becomes increasingly precise as J becomes large. Then 

[° pr(o) da = lim , (6) 

J - OD y— OO J 
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,V( * ) is the number of cases for which all elements of o are less 
than or equal to zero, that is. for which a m2X < 0. where is the 
maximum real eigenvalue component in o. An important feature of 
this definition is that it does not depend on the eigenvalues and 
eigenvectors retaining fixed structures. As parameters change, com- 
plex roots may coalesce to become real roots (or the reverse), and 
modes may exchange relative frequencies. The only matter for 
concern is whether or not all real pans of the eigenvalues remain in 
the left-half 5 plane. For J < », the probability of instability 
resulting from Monte Carlo evaluation is an estimate, denoted P. 

Histograms and cumulative distributions for varying degrees of 
stability are readily given by the Monte Carlo estimate of 
/f . pr( cr) do\ where 2 represents a maximum real eigenvalue 
component, and — oo < £ < oo. This histogram is a plot of N[(2 
- A) < < Z]/J versus 2; A is an increment in 2, jV[*] is 

the number of cases whose maximum real eigenvalue components 
lie in the increment, and J is the total number of evaluations. The 
histogram estimates the stability probability density function , 
pr(2), which is obtained in the limit for a continuous distribution of 
2 as A — 0 and J -* oo. The cumulative probability distribution 
of stability , Pr(2), is similarly estimated and presented as 
<2 )/ J versus 2, the exact distribution being achieved in the limit 
as J — ■ . Consequently 

? = 1 - Pr(0). (7) 

There is, of course, no limitation on admissible specifications for 
the multivariate pr(p); it may be Gaussian or non-Gaussian, as 
appropriate. Rayleigh, correlated, and any other well-posed distri- 
butions are admissible, the principal challenge being to properly 
shape (and correlate) the outputs of the random-number generator. 
In practice, system parameter uncertainties are most likely to be 
bounded, as typical quality control procedures eliminate out-of- 
tolerance devices, and there are physical limitations on component 
size, weight, shape, etc. The uniform distribution is particularly 
interesting, as it readily models bounded uncertainty, and it is the 
default distribution of most algorithms for random-number genera- 
tion. Given binary distributions for each parameter, in which the 
elements of p take maximum or minimum values with equal 
probability, the Monte Carlo evaluation reduces to 2 r deterministic 
evaluations, the result is exact, and the probability associated with 
each possible value of p is 1/2C Similarly, the distribution for r 
parameters, each of which takes w values (i.e., for quantized 
uniform distributions), can be obtained from w r evaluations; the 
probability of acquiring each value of p is 1/ w r . 

When has stochastic robustness been achieved? The answer is 
problem-dependent. pr(p) should be chosen to reflect physical 
limits of parameter uncertainty. In some applications involving 
bounded parameters, it will be possible to choose C such that 

= 0, and that is a desirable goal; however, if admissible parame- 
ter variations are unbounded, if C is constrained, or if the rank of 
CH is less than n, the minimum ~P may be greater than zero. C 
then must be chosen to satisfy performance goals and one of two 
robustness criteria; minimum ?, or P small enough to meet a 
reliability specification (e.g., one chance of instability in some 
larger number of realizations). One may object that the parameter 
distributions must be known or estimated for stochastic robustness 
analysis. However, if robustness estimates are strongly dependent 
on the statistics, then it is incumbent on the designer to know- 
something about the statistics. Otherwise, the final control system 
may be either unduly conservative (at the expense of performance) 
or insufficiently robust in the face of reai-worid uncertainties. 

II. Stochastic Root Locus 

While it is not necessary’ to plot the eigenvalues of (4) to 
determine or portray stochastic robustness, stochastic root loci 
provide insight regarding the effects of parameter uncertainty on 
system stability. Consider, for example, a classical second-order 
system with characteristic equation 

S 1 -I- 2f U n s + «; = 0. (8) 



(a) 



ing ratio and natural frequency. = 0.707, = 1; 50 CKK) Monte Carlo 

evaluations. The root density pr (X) is given in units of roots /unit length 
along the real axis and roots /unit area in the complex plane, (a) Scatter plot, 
(b) Oblique three-dimensional representation. 


Suppose that the damping ratio (f) and natural frequency («„) are 
nominally 0.707 and 1, respectively, and that each may be a 
Gaussian-distributed random variable with standard deviation of 
0.2. Root loci for individual parameter variations would follow 
classical configurations of root locus construction, with the heaviest 
density of roots in the vicinities of the nominal roots. If both f and 
Q)„ are uncertain and uncorrelated (i.e., p = the possible 

root locations become “clouds” surrounding the nominal values 
[Fig. 1(a)]. Understanding of robustness issues can be gained by 
plotting the density of the roots in a third dimension above the root 
locus plot [Fig. 1(b)]. This can be done by simply dividing the s 
plane into subspaces (or “bins”) and counting the number of roots 
in each bin as a sampled estimate of the root density p. The result 
is a multivariate histogram, with a and u> serving as independent 
variables. Complex root bins are elemental areas, for which p A is 
defined in units of roots /unit area. Real root bins are confined to the 
real axis; hence, p L measures roots /unit length. The density of 
roots depicts the likelihood that eigenvalues vary from their nominal 
values, including branches on the real axis and in the right- half s 
plane for large enough variations of f and w„. Fig. 1 is based on 
50 (XX) Monte Carlo evaluations, and numerical smoothing has been 
applied to account for sampling effects. An example of the his- 
togram and cumulative distribution is given in Fig. 2. The probabil- 
ity -of- instability estimate (2 = 0 on the cumulative distribution) is 
8 x 10“\ 

When considering instability, distinction must be made between 
the number of cases with right-half plane roots and the number of 
roots in the right-half plane. For example, a third-order system can 
be unstable » ith 1, 2, or 3 roots in the right-half plane, yet N 
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Fig. 2. (a) Histogram and (b) cumulative distribution for the second-order 

system with Gaussian damping rauo and natural frequency. 

would be incremented by one in each case. A high-order system 
with real roots could be unstable with one or more roots in the same 
right-half plane bin. Again, N would be incremented by one, 
although the bin's p depends upon the number of roots it contains. 

m. Confidence Intervals 

Because it is an estimate. P must have associated with it some 
notion of accuracy or relationship to the true^ underlying probability 
of instability P. Confidence intervals relate P to P by bounding the 
(unknown) true value with defined certainty. A confidence statement 
for P is 

Pr(I < P < U) - 1 - a (9) 

where ( L, U) is the interval estimate (lower and upper bounds), 

1 - cs is the confidence coefficient, and P lies within (L, U) with 
100(1 - cl)% confidence. Appropriate selection of the confidence 
coefficient and the number of evaluations determines the interval 
width. 

The method used to compute confidence intervals depends on the 
underlying probability distribution of the variable in question. The 
probability of instability is a binomial variable , with the outcome 
of a trial taking one of two possible values (stable or unstable) for 
each Monte Carlo evaluation. P therefore has a binomial probability 
density and cumulative distribution given by 

PK*) = - S>) J ~ X (10) 

x) = Y. (11) 

j^X 

where x is the number of occurrences of instability in J evalua- 
tions, and (/, x) is the binomial coefficient, J\( x\{J - x)l. The 
binomial test is applied to determine exact confidence intervals for 
binomial variables. The lower and upper confidence bounds are 
derived from the cumulative distribution and satisfy [19] 

Pt(X<x - 1 ) = £ - L) J ~ J = 1 - 2 ( 12 ) 

j-0 2 

Pr(X£x) = £ (J.J)C /'(1 - = ?■ (13) 

7*0 2 

where 1 — a is the required confidence coefficient and AMs the 
actual number of unstable cases after J evaluations ( X = jf>). 

The validity of the Monte Carlo analysis depends on a number of 
simulation parameters: the number of eigenvalues computed, the 
number of varying parameters, and their probability distributions. 
However, by applying the binomial test, the derivation of explicit 


rj 

! 


Pr(I) 


J-il. 



-t.O 0.)) 1.0 2.0 


I 

(b) 


I 

£ 

Z 


lc+05 


te*03 


le-Ol 



lc-06 le-05 le-04 0.00! 0.01 0.1 : 1 

0.5 

P (0 < P £ 0.5). or II * PI {0.3 < Pi 1) 


Fig. 3. Number of evaluations required for given confidence interval 
widths and confidence coefficient 1 - a = 0.95. Interval width is given as 
percent of p. For P > 0.5. the curves are symmetric, with the number of 
evaluations given by that for 1 - 


relationships between simulation parameters and the required num- 
ber of evaluations is avoided. Nevertheless, it offers a rigorous 
theory by which to calculate exact confidence intervals. The number 
of evaluations required for a specified confidence interval can be 
related to a single variable, P, and this relationship is valid for any 
simulation parameters or application. 

Fig. 3 presents the number of evaluations required for specified 
interval widths and a 95% confidence coefficient, given as a percent- 
age of P. An estimate P based on a small number of evaluations can 
be used as the abscissa of Fig. 3 to forecast the total number of 
evaluations required for the desired interval width. A larger confi- 
dence coefficient shifts the curves of Fig. 3 to higher numbers of 
evaluations. 


IV. Stochastic Robustness Example 

An example of the application of stochastic robustness is based on 
the longitudinal dynamics and control of an open- loop-unstable 
aircraft. The Forward-Swept- Wing Demonstrator’s aerodynamic 
center is forward of its center of gravity, resulting in static instabil- 
ity. Its stability matrix, control-effect matrix, and open-loop eigen- 
values are 


F = 


- 0.02 

- 0.0001 

0 

0 


G — 


\.4 = -0.1 


-0.3 

-0.4 

- 32.2 * 


- 1.2 
18. 

1 

-0.6 

0 

0 

(14) 

0 

I 

0 


~ - .04 

35. 1 



0 

0.2 

0 

-0.2 


(15) 

0 

0 



± 0.057/, 

- 5.15, 3.35. 

(16) 


The state components represent forward velocity, angle of attack, 
pitch rate, and pitch angle. The principal control surfaces are the 
canard control surface and the thrust setting. Possible uncertainties 
in aerodynamic and thrust effects as well as separate dynamic 
pressure (p and V) effects lead to a 12-element parameter vector 
(the remaining terms are kinematic, due to gravity, identically zero, 
or otherwise negligible) 


P “ [ P ^/n/ 12 / 13 / 22/32 &12&31 £ 32 ] * (l?) 

Although velocity (V) and air-density (p) are essentially determinis- 
tic, including them as parameters gives the ability to look at flight 
condition perturbations around the nominal value and reduces corre- 
lation of the remaining parameter. . p and V are modeled as uniform 
parameters, giving an indication of stochastic robustness over a 
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Case (a) 

C = | 

'0.1714 
. 0.984 

130.26 

-11.387 

33.165 

-2.968 

-l'm] G-diagU. 1.1.0) = diag (1,1) X = - 

35.0, - 5.14, 

- 3.32, - 

0.0183 

Case (b) 

I 

c = 

f 0.0270 
L 0.0107 

82.659 

-62.623 

20.927 
- 16.203 

— ° 902 8 ] 2 = dia S(1.1.1.0> R = 1000 diag (1, 1) 

X = -5.15, - 

- 3.36, - 

1.09, - 0.0186 

Case (c) 

c = 

f 0.1349 
l 0.0535 

413.294 

-313.112 

104.633 

-81.015 

~? ^ 1 1 X = -32.21, - 5.15, - 3.44, - 0.0! 
-9,509 J 





range of flight conditions. In terms of the elements p, i\ and G are 




-2gf lx pV 2 f xl 


V 
— 45 
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fiVfzz 

2 

pV 2 U_ 

2 
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G = 


pV 2 

2 


pV n 


pVfy. 


Sn 

*12 

0 

0 

* 31 

*32 

0 

0 


(18) 


(19) 


Linear-quadratic controllers are designed according to three spec- 
ifications: (a) Q ss diag (1, 1, 1,0), R — diag (1, 1); (b) Q = diag 


(1, 1, 1,0), M = diag (10CK), 1(XX)); and (c) the control gain matrix 
of Case (b) is multiplied by an arbitrary factor (5) to restore the 
closed-loop bandwidth to that of Case (a). The resulting control gain 
matrices and corresponding nominal closed-loop eigenvalues are 
given in Table I. These three cases have not been chosen to satisfy 
any particular performance criteria, but merely to demonstrate the 
impact of differing generalized design criteria on stochastic robust- 
ness. Furthermore, the designs are not meant to reflect acceptable 
control laws, as the high gams were purposely chosen to magnify 
robustness problems and to illustrate the application of stochastic 
robustness. 

For illustration, p and V are ±30% uniform parameters, and the 
remaining elements of p are subject to independent 30% standard 
deviation Gaussian uncertainties. Fig. 4 shows the stochastic root 
loci for each case, based on 25 (XX) Monte Carlo evaluations. In 
each case, the eigenvalue near the origin is least affected by the 
parameter changes, and its peak dominates the distribution. In Cases 
(a) and (c), the eigenvalue farthest to the left (not shown in Fig. 4) 
has an enc nmous variance along the real axis. Interaction »f roots 
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Fig. 5. 95% confidence intervals (a = 0.05) based on the binomial test for 
the Forward-Swept- Wing Demonstrator Aircraft with Gaussian parameters. 
Cases (a), (b), and (c). Solid line represents ? estimate, and dashed lines 
give confidence intervals. 

around the origin causes instability. The corresponding probability- 
of-instability estimates and 95% confidence intervals are (a) 0.0724 
(0.0692, 0.0756), (b) 0.0205 (0.0187,0.0222), and (c) 0.0076 
(0.0Q65, 0.0086). Robustness improves from Case (a) to (b) as 
control usage is restrained by high control weighting, and the ad hoc 
robustness recovery technique used in Case (c) gives additional 
improvement. Fig. 5 shows 95% confidence intervals and the P 
history with the number of evaluations. While confidence intervals 
for Cases (b) and (c) initially overlap, 25 (XX) evaluations are more 
than sufficient to rank the three cases in order of robustness. 

For uniformly distributed parameters in [0.7p, 1.3p], the extent 
of the parameter and eigenvalue distributions decreases substantially 
(Fig. 6). A comparison of Figs. 4 and 6 gives a better indication of 
the effects of Gaussian “tails” on the eigenvalue probability densi- 
ties. J> and confidence intervals for 25 (XX) evaluations are (a) 
3.6E-4 (1.25E^, 5.95E^) and (b) and (c) zero (0, 1.48E-4). For 
12-parameter binary distributions, 4096 evaluations produce exact 
results. ? for each case is (a) 0.1191 and (b) and (c) zero. 

Stochastic stability robustness is given as a function of the control 
design parameter v, where the LQR control weighting matrix 
R = vl 2 (Fig. 7). Under the specified limits of parameter uncer- 
tainty, the distinction in stability and robustness versus design 
parameter statistics is apparent. Note that the qualitative result that 
10 < v < 10 5 provides the most robust designs is independent of 
the assumed probability distribution. In particular, the fact that the 
minimum-control-energy case {y — oo) represents the least robust 
design would not be obvious using standard robustness analysis 
techniques such as unstructured-singular-value analysis. While the 
exact relationship between parameter uncertainty and eigenvalue 
location is complicated, the increase in P beyond 10 5 is attributed to 
the shift in the closed-loop eigenvalues towards the imaginary axis. 
For example, at v = I0 8 the closed-loop eigenvalues are -5.15, — 
3.35, — 0.0104 ± 0.057y; the first two are almost identical to those 
of Case (b), while the complex pair is closer to the imaginary axis. 
The kind of results presented in Fig. 7 offer controller design insight 
and show nonobvious robustness characteristics. 

V. Conclusion 

Stochastic robustness offers a rigorous yet straightforward alterna- 
tive to current metrics for control system robustness that is simple to 
compute and is unfettered by normally difficult problem statements, 
such as non-Gaussian statistics, products of parameter variations, 
and structured uncertainty. The approach answers the question, 
“How likely is the closed-loop system to fail, given limits of 
parameter uncertainty? “ It makes good use of modem computa- 
tional and graphic tools, and it is easily related to practical design 
considerations. The principal difficulty in applying this method to 
controlled systems is that it is computationally intensive; however, 
requirements are well within the capabilities of existing computers. 
The principal advantage of the approach is that it is easily imple- 
mented, and results have direct bearing on engineering objectives. 


' pn'X) 



Aircraft with uniformly distributed parameters. Cases (a), (b), and (c). The 
root density pr (X) is given in units of roots /unit length along the real axis 
and roots /unit area in the complex plane.' 
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Fig. 7. Probability -of- instability estimates versus LQR design parameter v 

( R = v/ a ) for forward-swept wirg aircraft example. 
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